# Replication materials for 
# Peisker, Hoffmann, Muttarak (2025)
# Climate news mediates extreme weather effects on climate change concern

packages <- c(
  "tidyverse", "broom", "data.table", "fixest", "parallel",
  "texreg", "scico", "marginaleffects", "patchwork"
  )
lapply(packages, library, character.only = TRUE)
lapply(packages, citation)

source("R Files/functions.R")
theme_set(theme_bw())
setDTthreads(parallel::detectCores())
options("mc.cores" = detectCores())

#### load data and set variable names ####
load("Data/eb_climate_reg.RData")
eb_panel <- panel(eb_climate_weekly_n, ~unit+t)
eb_panel_outcome <- panel(eb_climate_outcome, ~unit+t)

#### set variable names ####
# independent variables
names(eb_climate_weekly_n)
temp_4w <- c(
  "tmean_nf_4w",
  "tmean_nf_neg1_4w + tmean_nf_pos1_4w",
  "tmean_hs_q05_4w + tmean_hs_q95_4w",
  "tmean_xf_q05_4w + tmean_xf_q95_4w"
)
emm_l1 <- c(
  "emm_climate_w_reg_l1", "emm_climate_w_nat_l1", 
  "emm_env_w_reg_l1", "emm_env_w_nat_l1" 
)

# controls
cpol_vars <- 
  grep("_lon", names(eb_climate_weekly_n), value = T) %>% 
  .[!grepl("2014", .)]
lagged_media <- " + emm_climate_w_nat_4w_l2 + emm_climate_w_reg_4w_l2 + emm_env_w_nat_4w_l2 + emm_env_w_reg_4w_l2 + emm_energy_w_nat_4w_l2 + emm_energy_w_reg_4w_l2"

# dictionary of variable names
temp_4w_map <- list(
  "tmean_nf_4w" = "Temperature anomaly", 
  "tmean_nf_pos1_4w" = "Positive temp. anomaly (1 SD)", 
  "tmean_nf_neg1_4w" = "Negative temp. anomaly (1 SD)", 
  "tmean_hs_q975_4w" = "Warm spell severity (97.5\\ %)",
  "tmean_hs_q025_4w" = "Cold spell severity (2.5\\ %)", 
  "tmean_hs_q95_4w" = "Warm spell severity (95\\ %)",
  "tmean_hs_q05_4w" = "Cold spell severity (5\\ %)", 
  "tmean_hs_q90_4w" = "Warm spell severity (90\\ %)",
  "tmean_hs_q10_4w" = "Cold spell severity (10\\ %)",
  "tmean_xf_q95_4w" = "Warm spell duration (95\\ %)",
  "tmean_xf_q05_4w" = "Cold spell duration (5\\ %)", 
  "event_type_cop_l1" = "COP",  
  "event_type_sports_l1" = "Sports event", 
  "event_type_activism_l1" = "Activism",  
  "event_type_g7_l1" = "G7 summit",
  "event_type_year_l1" = "Year change",
  "event_type_eu_election_l1" = "EU election",
  "event_type_nat_election_l1" = "National election",
  "seasonspring" = "Spring",
  "seasonsummer" = "Summer",
  "seasonwinter" = "Winter",
  "tmean_nf_pos1_4w:event_type_cop_l1" = "Pos. temp. anomaly × COP", 
  "tmean_nf_pos1_4w:event_type_sports_l1" = "Pos. temp. anomaly × sports event",
  "tmean_nf_pos1_4w:seasonspring" = "Pos. temp. anomaly × spring", "tmean_nf_pos1_4w:seasonsummer" = "Pos. temp. anomaly × summer", "tmean_nf_pos1_4w:seasonwinter" = "Pos. temp. anomaly × winter",
  "tmean_nf_neg1_4w:event_type_cop_l1" = "Neg. temp. anomaly × COP", 
  "tmean_nf_neg1_4w:event_type_sports_l1" = "Neg. temp. anomaly × sports event",
  "tmean_nf_neg1_4w:seasonspring" = "Neg. temp. anomaly × spring", "tmean_nf_neg1_4w:seasonsummer" = "Neg. temp. anomaly × summer", "tmean_nf_neg1_4w:seasonwinter" = "Neg. temp. anomaly × winter",  "emm_climate_w_reg_l1" = "Regional climate news", 
  "emm_climate_w_reg_l1" = "Regional climate news",
  "emm_climate_w_nat_l1" = "National climate news",
  "emm_env_w_reg_l1" = "Regional environmental news", 
  "emm_env_w_nat_l1" = "National environmental news",
  "emm_energy_w_reg_l1" = "Regional energy news", 
  "emm_energy_w_nat_l1" = "National energy news",
  "emm_climate_w_reg_4w_l2" = "Regional climate news",
  "emm_climate_w_nat_4w_l2" = "National climate news",
  "emm_env_w_reg_4w_l2" = "Regional environmental news",
  "emm_env_w_nat_4w_l2" = "National environmental news",
  "emm_energy_w_reg_4w_l2" = "Regional energy news", 
  "emm_energy_w_nat_4w_l2" = "National energy news"
  )
coef_map_outcome <- list(
  "tmean_nf_pos1_4w" = "$T$: Pos. temperature anomaly",
  "tmean_nf_neg1_4w" = "Neg. temperature anomaly", 
  "emm_climate_w_nat_4w" = "$M$: National climate news", 
  "emm_climate_w_reg_4w" = "Regional climate news", 
  "emm_env_w_nat_4w" = "National environmental news", 
  "emm_env_w_reg_4w" = "Regional environmental news", 
  "emm_energy_w_nat_4w" = "National energy news", 
  "emm_energy_w_reg_4w" = "Regional energy news", 
  "google_trend_4w" = "Google web searches"
)
coef_map_outcome_poly <- list(
  "tmean_nf_pos1_4w" = "Pos. temperature anomaly",
  "tmean_nf_neg1_4w" = "Neg. temperature anomaly", 
  "emm_climate_w_nat_4w" = "National climate news", 
  "I(emm_climate_w_nat_4w^2)" = "National climate news²",
  "I(emm_climate_w_nat_4w^3)" = "National climate news³",
  "emm_climate_w_reg_4w" = "Regional climate news", 
  "I(emm_climate_w_reg_4w^2)" = "Regional climate news²",
  "I(emm_climate_w_reg_4w^3)" = "Regional climate news³",
  "emm_env_w_nat_4w" = "National environmental news", 
  "emm_env_w_reg_4w" = "Regional environmental news", 
  "emm_energy_w_nat_4w" = "National energy news", 
  "emm_energy_w_reg_4w" = "Regional energy news", 
  "google_trend_4w" = "Google web searches"
)
gof_baseline = c(
  "Observations", "Region-year fixed effects (groups)", "Season fixed effects (groups)", 
  "Event type fixed effects (groups)", "R² (overall)", "R² (within)"
)

#### mediator equations, outcome: media coverage ####
#### baseline models ####
yv <- emm_l1
xv <- temp_4w

tab_emm_baseline <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_baseline)))
tab_export <- tab_emm_baseline[col_export]
screenreg(tab_export)
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_baseline.tex",
  custom.header = list("Regional climate news" = 1:4, "National climate news" = 5:8),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
  )

#### dot-whisker baseline ####
tab_emm_baseline_tidy <- 
  lapply(tab_emm_baseline, function(x){cbind(lhs = as.character(x$fml[[2]]), tidy(x))}) %>% 
  bind_rows() %>% 
  mutate(
    ci_low = estimate - 1.96 * std.error,
    ci_high = estimate + 1.96 * std.error,
    Topic = case_when(
      grepl("climate", lhs) ~ "Topic: Climate change",
      grepl("env", lhs) ~ "Topic: Environment & Ecology"
    ),
    Publication = factor(case_when(
      grepl("reg", lhs) ~ "Regional climate news",
      grepl("nat", lhs) ~ "National climate news"
    ), levels = c("Regional climate news", "National climate news")
    ),
    Type = if_else(grepl("(q05)|(_neg1)", term), "Cold spell", "Warm spell"),
    var = factor(case_when(
      grepl("hs", term) ~ "Spell severity",
      grepl("xf", term) ~ "Spell duration",
      grepl("nf_4w", term) ~ "Temperature anomaly",
      grepl("1_4w", term) ~ "Pos./Neg. temperature anomaly"
    ), levels = rev(c(
      "Temperature anomaly", "Pos./Neg. temperature anomaly", "Spell severity", "Spell duration"
    )))
  )

tab_emm_baseline_tidy %>%
  filter(Topic == "Topic: Climate change", !is.na(var)) %>% 
  ggplot(aes(y = var, x = estimate)) +
  geom_vline(xintercept = 0, colour = "grey60") +  
  geom_pointrange(
    aes(xmin = ci_low, xmax = ci_high, shape = Type),
    position = position_dodge(width = 0.5)
  ) +
  facet_wrap(vars(Publication)) +
  scale_color_scico_d(palette = "vik", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = "Standardized effect on climate news") + 
  theme(
    legend.position = "bottom", 
    legend.title = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_text(angle = 45, vjust = 0.5, hjust=1)
    ) 
ggsave("Figures/5dw_emm_baseline.pdf", width = 8, height = 4)

#### controls for spatio-temporal polynomials ####
yv <- emm_l1
xv <- paste0(temp_4w, " + ", paste0(cpol_vars, collapse = " + "))

tab_emm_sp <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = c("unit", "year", "season", "event")),#
    eb_panel, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_sp)))
tab_export <- tab_emm_sp[col_export]
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_spat_poly.tex",
  custom.header = list("Regional climate news" = 1:4, "National climate news" = 5:8),
  custom.coef.map	= temp_4w_map,
  custom.gof.row = list("Spatio-temporal polynomials" = rep("Yes", 8)),
  custom.gof.names = c("Observations", "Region fixed effects (groups)", "Year fixed effects (groups)", 
                       "Season fixed effects (groups)", "Event type fixed effects (groups)", "R² (overall)", "R² (within)")
)

#### with conley standard errors ####
yv <- emm_l1
xv <- temp_4w

tab_emm_conley <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel, vcov = conley(500)
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_conley)))
tab_export <- tab_emm_conley[col_export]
my_screenreg(tab_export, custom.coef.map = temp_4w_map)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_conley.tex",
  custom.header = list("Regional climate news" = 1:4, "National climate news" = 5:8),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
)

#### unweighted articles ####
yv <- c("emm_climate_reg_l1", "emm_climate_nat_l1")
xv <- temp_4w

tab_emm_w <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_climate_weekly_n, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_w)))
tab_export <- tab_emm_w[col_export]
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_w.tex",
  custom.header = list("Regional climate news" = 1:4, "National climate news" = 5:8),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
)

#### different thresholds ####
xv <- c(
  "tmean_hs_q975_4w + tmean_hs_q025_4w",
  "tmean_hs_q95_4w + tmean_hs_q05_4w",
  "tmean_hs_q90_4w + tmean_hs_q10_4w"
)
yv <- emm_l1

tab_emm_thr <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_climate_weekly_n, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_thr)))
tab_export <- tab_emm_thr[col_export]
my_screenreg(tab_export, custom.coef.map	= temp_4w_map)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_thr.tex",
  custom.header = list("Regional climate news" = 1:3, "National climate news" = 4:6),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
)

#### controls for lagged media ####
yv <- emm_l1
xv <- paste0(temp_4w, lagged_media)

tab_emm_lag <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel, vcov = cluster ~ nuts0year 
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

col_export <- which(grepl("climate", names(tab_emm_lag)))
tab_export <- tab_emm_lag[col_export]
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_emm_lag.tex",
  custom.header = list("Regional climate news" = 1:4, "National climate news" = 5:8),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
)

#### distributed lag models ####
emm_ind_vars <- c(
  "Warm spell" = "tmean_nf_pos1",
  "Cold spell" = "tmean_nf_neg1", 
  "COP" = "event_type_cop",
  "G7 summit" = "event_type_g7",
  "EU election" = "event_type_eu_election",
  "National election" = "event_type_nat_election",
  "Sports event" = "event_type_sports"
)
emm_dep_vars <- c(
  "Regional outlets" = "_reg",
  "National outlets" = "_nat"
)
lags <- -2:8
tab_emm_dl <- feols(
  c(emm_climate_w_nat, emm_climate_w_reg) ~ 
    l(tmean_nf_pos1, lags) + l(tmean_nf_neg1, lags) + 
    l(event_type_cop, lags) + l(event_type_g7, lags) + 
    l(event_type_sports, lags) + 
    l(event_type_eu_election, lags) + l(event_type_nat_election, lags) + 
    l(emm_climate_w_nat, 1:8) + l(emm_climate_w_reg, 1:8) +
    l(emm_env_w_nat, 1:8) + l(emm_env_w_reg, 1:8) +
    l(emm_energy_w_nat, 1:8) + l(emm_energy_w_reg, 1:8) +
    event_type_year  
  | nutsyear + season, # FE
  vcov = cluster ~ nuts0year,
  data = eb_panel
)
make_lc_plot(tab_emm_dl) %>%
  plot_lc("Figures/5emm_dist_lags.pdf")


#### outcome: google trend ####
yv <- c("google_trend_l1")
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "emm_climate_w_reg_l1 + emm_climate_w_nat_l1 + emm_env_w_reg_l1 + emm_env_w_nat_l1 + emm_energy_w_reg_l1 + emm_energy_w_nat_l1",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_l1 + emm_climate_w_nat_l1 + emm_env_w_reg_l1 + emm_env_w_nat_l1 + emm_energy_w_reg_l1 + emm_energy_w_nat_l1"
)

tab_google_emm <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel, vcov = cluster ~ nuts0year #conley(500)# 
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- tab_google_emm
my_screenreg(tab_export, custom.coef.map	= temp_4w_map)
my_texreg(
  tab_export,
  file = "Tables/2tab_google_emm.tex",
  custom.header = list("Google web searches for climate change" = 1:length(tab_export)),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline
)

#### placebo tests mediator ####
set.seed(123)
shuffle_vars <- c("tmean_nf_neg1_4w", "tmean_nf_pos1_4w")#
out_shuffle <- paste0(shuffle_vars, "_shuffle")
sel_vars <- c(shuffle_vars, m_fv, "emm_climate_w_nat_l1", "nuts0year", "lat", "lon")
placebo_mediator <- replicate(
  1e4,
  feols(
      f_feols(yvar = "emm_climate_w_nat_l1", xvars = out_shuffle, fevars = m_fv), 
      data = (eb_climate_weekly_n[, ..sel_vars, with = FALSE] %>% 
        .[,
          c(out_shuffle) := lapply(.SD, sample),
          by = .(nutsyear),
          .SDcols = c(shuffle_vars)
        ]),
      vcov = cluster ~ nuts0year
    ) %>% 
    .$coeftable,
  simplify = FALSE
)
# save(placebo_mediator, file = "Data/placebo_mediator.RData")

# load("Data/placebo_mediator.RData")
placebo_mediator_data <- 
  placebo_mediator %>%
  lapply(as.data.frame) %>% 
  bind_rows %>%
  rownames_to_column(var = "var") %>% 
  mutate(
    var = case_when(
      grepl("tmean_nf_pos1_4w", var) ~ "Positive temperature anomaly",
      grepl("tmean_nf_neg1_4w", var) ~ "Negative temperature anomaly",
      )
  ) 

vline_df <-
  tab_emm_baseline[which(grepl("climate", names(tab_emm_baseline), fixed = TRUE))] %>% 
  .[which(grepl("nat", names(.), fixed = TRUE))] %>% 
  .[[2]] %>% .$coeftable %>% as_tibble %>% 
  cbind(var = c("Negative temperature anomaly", "Positive temperature anomaly"))
sd_df <-  
  placebo_mediator_data %>% 
  group_by(var) %>% 
  summarise(sd = sd(Estimate))
placebo_coef <- 
  placebo_data %>% 
  left_join(sd_df) %>% 
  mutate(
    norm_dist = dnorm(Estimate, mean = 0, sd = sd)
  ) %>% 
  ggplot() +
  geom_vline(xintercept = 0) +
  geom_vline(data = vline_df, aes(xintercept = Estimate), color = "red") +
  geom_vline(data = vline_df, aes(xintercept = Estimate + 1.96 * `Std. Error`), color = "red", linetype = "dashed") +
  geom_vline(data = vline_df, aes(xintercept = Estimate - 1.96 * `Std. Error`), color = "red", linetype = "dashed") +
  geom_line(aes(x = Estimate, y = norm_dist), color = "blue") +
  geom_density(aes(Estimate), adjust = 2, trim = TRUE) +
  facet_wrap(~var) +
  labs(y = "density")
# placebo_coef

placebo_t <- 
  placebo_data %>% 
  ggplot() +
  geom_vline(xintercept = 0) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1), color = "blue") +
  geom_density(aes(`t value`), adjust = 2, trim = TRUE) +
  facet_wrap(~var)+
  labs(y = "density")
placebo_coef / placebo_t
ggsave("Figures/5placebo_mediator.pdf", width = 7, height = 8)

#### interactions with season, ipcc, and sports ####
yv <- emm_l1
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + 
  event_type_sports_l1 + event_type_cop_l1 + event_type_eu_election_l1 + event_type_nat_election_l1 + 
  event_type_g7_l1 + event_type_year_l1 + event_type_activism_l1 + season +
  tmean_nf_pos1_4w:event_type_cop_l1 + tmean_nf_pos1_4w:event_type_sports_l1 + tmean_nf_pos1_4w:season +
  tmean_nf_neg1_4w:event_type_cop_l1 + tmean_nf_neg1_4w:event_type_sports_l1 + tmean_nf_neg1_4w:season"
)

tab_emm_int <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = c("nutsyear")),
    eb_panel, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- tab_emm_int
my_screenreg(tab_export)

my_texreg(
  tab_export,
  file = "Tables/2tab_emm_int.tex",
  custom.header = list("Regional climate news" = 1, "National climate news" = 2, "Regional environmental news" = 3, "National environmental news" = 4),
  custom.coef.map	= temp_4w_map,
  custom.gof.names = gof_baseline[!grepl("Event|Season", gof_baseline)]
)

#### calc linear combinations for COP and sports ####
coef(tab_emm_int[[2]]) %>% names %>% tibble() %>% print(n=50)
hypo_contr <- matrix(
  c(1, 0, rep(0, 10), 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
    1, 0, rep(0, 10), 1, 0, 0, 0, 1, 0, 0, 0, 0, 0,
    1, 0, rep(0, 10), 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
    1, 0, rep(0, 10), 0, 1, 0, 1, 0, 0, 0, 0, 0, 0
    ),
  byrow = FALSE, nrow = coef(tab_emm_int[[1]]) %>% length,
  dimnames = list(
    tab_emm_int[[1]] %>% coef %>% names,
    c("no_cop", "cop", "no_sports", "sports"))
  )
hypo_contr

tab_lincom <- 
    cbind(
      lhs = tab_emm_int[[1]]$fml[[2]] %>% as.character,
      hypotheses(tab_emm_int[[1]], hypo_contr)
    ) %>% 
  bind_rows(
    cbind(
      lhs = tab_emm_int[[2]]$fml[[2]] %>% as.character,
      hypotheses(tab_emm_int[[2]], hypo_contr)
    ) 
  ) %>% 
  mutate(
    event = factor(
      if_else(grepl("no", term), "Without event", "With event"),
      levels = c("Without event", "With event")),
    term = factor(case_when(
      grepl("cop", term) ~ "COP (winter)",
      grepl("sports", term) ~ "Sports (summer)"
    ), levels = c("Sports (summer)", "COP (winter)")),
    publication = factor(case_when(
      lhs == "emm_climate_w_reg_l1" ~ "Regional climate news",
      lhs == "emm_climate_w_nat_l1" ~ "National climate news"
      ), levels = c("Regional climate news", "National climate news"))
    )
tab_lincom

tab_lincom %>% 
  ggplot(aes(y = term, x = estimate)) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high, shape = event),
    position = position_dodge(width = -0.75)
  ) +
  geom_vline(xintercept = 0, colour = "grey60") +  
  facet_wrap(vars(publication)) +
  scale_color_scico_d(palette = "vik", direction = 1, begin = 0.2, end = 0.8) +
  labs(x = "Standardized effect on climate news") + 
  theme(
    legend.position = "bottom", 
    legend.title = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_text(angle = 45, vjust = 0.5, hjust=1)
  ) 
ggsave("Figures/5dw_emm_lincom.pdf", width = 7, height = 3)

#### Table 1 ####
yv <- "emm_climate_w_nat_4w"
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
)
tab_mediator <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel_outcome[!is.na(cc_eu_concern_w),], vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)
screenreg(tab_mediator)

yv <- "cc_eu_concern_w"
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
)
tab_outcome <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel_outcome, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- c(tab_mediator, tab_outcome)
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_mediate.tex",
  custom.header = list("$M$: National climate news" = 1:3, "$Y$: Climate change concern" = 4:7),
  custom.coef.map	= coef_map_outcome, 
  custom.gof.names = gof_baseline
)

#### control for spatio-temporal trends ####
yv <- "emm_climate_w_nat_4w"
xv <- paste0(c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
), " + ", paste0(cpol_vars, collapse = " + "))

tab_mediator_sp <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = c("unit", "year", "season", "event")), 
    eb_panel_outcome[!is.na(cc_eu_concern_w),], vcov = cluster ~ nuts0year 
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)
screenreg(tab_mediator_sp)

yv <- "cc_eu_concern_w"
xv <- paste0(c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
), " + ", paste0(cpol_vars, collapse = " + "))
tab_outcome_sp <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = c("unit", "year", "season", "event")), 
    eb_panel_outcome, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- c(tab_mediator_sp, tab_outcome_sp)
my_screenreg(tab_export, omit.coef = "lat")
my_texreg(
  tab_export,
  file = "Tables/2tab_mediate_sp.tex",
  custom.header = list("$M$: National climate news" = 1:3, "$Y$: Climate change concern" = 4:7),
  custom.coef.map	= coef_map_outcome,
  custom.gof.row = list("Spatio-temporal polynomials" = rep("Yes", length(tab_export))),
  custom.gof.names = c("Observations", "Region fixed effects (groups)", "Year fixed effects (groups)", 
                       "Season fixed effects (groups)", "Event type fixed effects (groups)", "R² (overall)", "R² (within)")
)

#### non-linear media effects ####
yv <- "cc_eu_concern_w"
xv <- c(
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_reg_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_nat_4w^2 + emm_climate_w_reg_4w + emm_climate_w_reg_4w^2",
  "tmean_nf_pos1_4w + tmean_nf_neg1_4w + emm_climate_w_nat_4w + emm_climate_w_nat_4w^2 + emm_climate_w_reg_4w + emm_climate_w_reg_4w^2 + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w + google_trend_4w"
)
tab_outcome_poly <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel_outcome, vcov = cluster ~ nuts0year
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- tab_outcome_poly
my_screenreg(tab_export, custom.coef.map	= coef_map_outcome_poly)
texreg(
  tab_export,
  file = "Tables/2tab_concern_poly.tex",
  custom.header = list("Climate change concern" = 1:length(tab_export)),
  custom.coef.map	= coef_map_outcome_poly, 
  # custom.gof.names = gof_baseline, 
  include.adjrs = TRUE, include.fstatistic = FALSE, include.groups = TRUE,
  digits = 3, table = FALSE, booktabs = TRUE, use.packages = FALSE, 
  dcolumn = TRUE, stars = c(0.01, 0.05, 0.1),
  custom.model.names = paste0("(", 1:length(tab_export), ")")
)

plot_q_nat <- plot_slopes(
  tab_export[[4]],
  variables = "emm_climate_w_nat_4w",
  condition = "emm_climate_w_nat_4w"
) + 
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_vline(xintercept = 0, alpha = 0.5) +
  scale_y_continuous(limits = c(-0.5, 0.8)) +
  labs(x = "National climate news")
plot_q_reg <- plot_slopes(
  tab_export[[4]],
  variables = "emm_climate_w_reg_4w",
  condition = "emm_climate_w_reg_4w"
) + 
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_vline(xintercept = 0, alpha = 0.5) +
  scale_y_continuous(limits = c(-0.5, 0.8)) +
  labs(x = "Regional climate news")
plot_q_nat + plot_q_reg + 
  plot_annotation(tag_levels = "a")
ggsave("Figures/5concern_media_q.pdf", width = 8, height = 4)


#### placebo tests outcome ####
set.seed(123)
shuffle_vars <- c("tmean_nf_pos1_4w", "tmean_nf_neg1_4w", "emm_climate_w_nat_4w")#
out_shuffle <- paste0(shuffle_vars, "_shuffle")
sel_vars <- c(shuffle_vars, m_fv, "cc_eu_concern_w", "nuts0year", "lat", "lon")
placebo_outcome <- replicate(
  1e4,
  feols(
    f_feols(yvar = "cc_eu_concern_w", xvars = out_shuffle, fevars = m_fv), 
    data = (eb_climate_outcome[, ..sel_vars, with = FALSE] %>% 
              .[,
                c(out_shuffle) := lapply(.SD, sample),
                by = .(nutsyear),
                .SDcols = c(shuffle_vars)
              ]),
    vcov = cluster ~ nuts0year
  ) %>% 
    .$coeftable,
  simplify = FALSE
)
# save(placebo_outcome, file = "Data/placebo_outcome.RData")

# load("Data/placebo_outcome.RData")
placebo_outcome_data <- 
  placebo_outcome %>%
  lapply(as.data.frame) %>% 
  bind_rows %>%
  rownames_to_column(var = "var") %>% 
  mutate(
    var = case_when(
      grepl("tmean_nf_pos1_4w", var) ~ "Positive temperature anomaly",
      grepl("tmean_nf_neg1_4w", var) ~ "Negative temperature anomaly",
      grepl("emm_climate_w_nat_4w", var) ~ "National climate news"
    )
  ) 

vline_df <-
  tab_outcome[[2]]$coeftable %>% 
  as_tibble %>% 
  cbind(var = c("Positive temperature anomaly", "Negative temperature anomaly", "National climate news"))
sd_df <-  
  placebo_outcome_data %>% 
  group_by(var) %>% 
  summarise(sd = sd(Estimate))
placebo_coef <- 
  placebo_data %>% 
  left_join(sd_df) %>% 
  mutate(
    norm_dist = dnorm(Estimate, mean = 0, sd = sd)
  ) %>% 
  ggplot() +
  geom_vline(xintercept = 0) +
  geom_vline(data = vline_df, aes(xintercept = Estimate), color = "red") +
  geom_vline(data = vline_df, aes(xintercept = Estimate + 1.96 * `Std. Error`), color = "red", linetype = "dashed") +
  geom_vline(data = vline_df, aes(xintercept = Estimate - 1.96 * `Std. Error`), color = "red", linetype = "dashed") +
  geom_line(aes(x = Estimate, y = norm_dist), color = "blue") +
  geom_density(aes(Estimate), adjust = 2, trim = TRUE) +
  facet_wrap(~var) +
  labs(y = "density")
# placebo_coef

placebo_t <- 
  placebo_data %>% 
  ggplot() +
  geom_vline(xintercept = 0) +
  stat_function(fun = dnorm, args = list(mean = 0, sd = 1), color = "blue") +
  geom_density(aes(`t value`), adjust = 2, trim = TRUE) +
  facet_wrap(~var)+
  labs(y = "density")
placebo_coef / placebo_t
ggsave("Figures/5placebo_outcome.pdf", width = 9, height = 8)

#### reverse causality ####
yv <- c(
  "emm_climate_w_reg_4w_L1", "emm_climate_w_nat_4w_L1"
)
xv <- c(
  "cc_eu_concern_w",
  "cc_eu_concern_w + emm_climate_w_reg_4w + emm_climate_w_nat_4w + emm_env_w_reg_4w + emm_env_w_nat_4w + emm_energy_w_nat_4w + emm_energy_w_reg_4w"
)

tab_concern_lead <- mapply(
  function(y, x){  feols(
    f_feols(yvar = y, xvars = x, fevars = m_fv), 
    eb_panel_outcome, vcov = cluster ~ nuts0year 
  )},
  rep(yv, each = length(xv)), 
  rep(xv, times = length(yv)), 
  SIMPLIFY = FALSE
)

tab_export <- tab_concern_lead
my_screenreg(tab_export)
my_texreg(
  tab_export,
  file = "Tables/2tab_concern_lead.tex",
  custom.header = list("Regional climate news" = 1:2, "National climate news" = 3:4),
  custom.coef.names	= c(
    "Climate change concern",
    "Climate news (lag, regional)", "Climate news (lag, national)", 
    "Environmental news (lag, regional)", "Environmental news (lag, national)",
    "Energy news (lag, regional)", "Energy news (lag, national)"
  ),
  custom.gof.names = gof_baseline
)

